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We examine the hydrodynamics of a granular gas using numerical simulation. We demonstrate 
| the appearance of shearing and clustering instabilities predicted by linear stability analysis, and 

show that their appearance is directly related to the inelasticity of collisions in the material. We 
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discuss the rate at which these instabilities arise and the manner in which clusters grow and merge. 
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One of the key differences between a granular material and a regular fluid is that the grains of the former lose energy 
with each collision, while the molecules of the latter do not. Even when the inelasticity of the collisions is small, it 
can give rise to dramatic effects, such as the Maxwell Demon effect [1] and, the topic of this paper, the phenomenon of 
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granular clustering. Experiments[2, 3] and molecular dynamics simulations [4] alike show that low-density collections 
of grains ("granular gases") in the absence of gravity do not become homogeneous with time, but instead form denser 
clusters of stationary particles surrounded by a lower density region of more energetic particles. A kinetic explanation 
for this behavior is that, when a particle enters a region of slightly higher density, it has more collisions, loses more 
energy, and so is less able to leave that region. This increases the local density and makes it more likely that additional 
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system. Goldhirsch and Zanetti[4] , for instance, describe clustering as the result of a hydrodynamic instability: a region 



\ particles are captured in the same way. 

We are interested in describing this clustering behavior using hydrodynamics. There is considerable work[5] deriving 
granular hydrodynamics from kinetic theory, focusing on analytical treatments of the long- wavelength behavior of the 
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of slightly higher density has more collisions, so more energy is lost and the region has a lower "temperature" [6]. Less 

S temperature results in less pressure, and this lower pressure region in turn attracts more mass from the surrounding 
higher-pressure regions. Their paper uses long-wavelength stability analysis to show that, in a system of hydrodynamic 
equations similar to Eq. 1 below, higher-density regions do indeed have lower pressure, fueling the instability. 
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, In this paper, we study the hydrodynamics of granular clustering in zero gravity, by using numerical simulation. Our 



motivation is to determine whether a coarse-grained description, in terms of local particle, momentum, and energy 
densities, can be used to treat characteristic behaviors of granular materials as a self-contained dynamical system [7]. 
We show that the instabilities predicted by linear analysis do arise in our simulations, and discuss how the onset of 
these instabilities depends on the inelasticity of collisions in the material. We also show the manner in which clusters 
develop. 

We begin with a number density field p, a flow velocity field u, and a "temperature" [6] field T. These are related 
by a standard set of hydrodynamic equations for granular materials, introduced by Haff[8]: 
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where repeated indices are summed over, and where P is the pressure, k is the bare thermal conductivity, and 
Vijki = v(dik$ji + 5a$kj + SijSki) is the isotropic bare viscosity tensor. These equations bear much in common with 
those for normal fluids[9]. The most important addition is that of the term — 7T, which accounts for the inelasticity 
of collisions; the parameter 7 is proportional to (1 — r 2 ), where r is the coefficient of restitution. Using kinetic theory 
results [8], the transport coefficients are chosen to depend on temperature and density: 



Typically, work in granular hydrodynamics is done in low-density regimes, where grains may be treated as point 
particles interacting via collisions. When simulating aggregation, however, one must take excluded volume into 
account. We do this by introducing a barrier in the pressure P(p) at some maximum (close-packed) density p . This 
is in addition to the usual hydrodynamic pressure pT. We choose in particular the simple quadratic form 



where U is a positive parameter, 9{x) is the unit step function, and po is the close-packed density. This method, 
which we introduced in an earlier papcr[10], is a simple way to model the incompressibility of the system at high 
densities [11]. 

We evaluate our equations in two dimensions using a finite-difference Runge-Kutta method, on a square lattice with 
periodic boundary conditions. (See footnote [12] for more details.) The lattice spacing is chosen to be large enough 
so that each site contains a number of grains, and we can consider the density to be a continuous variable. We start 
with random initial conditions p = 0.1 + 0.001ri(z,x), u z — r2(z,x), u x = r^{z 1 x) 1 and T = 1 + O.lr^z, x), where 
ri{z,x) are random numbers chosen between —1 and 1. The model's other parameters for the data presented here 
are r]o = 25, kq = 1, U = 4 x 10 4 , /?o = 0.2, and 7 taking on several different values. All numbers given here are in 
dimensionless units[13]. Our time step in these units is At = 10~ 3 . 

We begin with a system that is 64 x 64 in size. The homogeneous state with which we initialize our system is 
already a solution to the equations above. In this initial homogeneous cooling state, the velocity and all gradients 
vanish, and the temperature decays with time due to the inelasticity. Eq.l reduces to 
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P = P T + U(p 2 - pl)6{p - po) 
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FIG. 1: The first plot shows the evolution of the average temperature of the system over time, for three different values of the 
inelasticity parameter 7, on a log- log plot. The second is the logarithmic derivative of T(t); that is, , or the slope of the 

lines in the first graph. The lower graph shows the exponent of the power-law decay rate. 

which yields Haff's cooling law, T(t) = T(0)(l + t/t )~ 2 - This state is seen universally in simulations[14-17], but only 
initially, for it is unstable to hydrodynamic modes [4], resulting in a long-range shear flow followed by the clustering 
instability mentioned above. 
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FIG. 2: The maximum power-law decay rate of the average temperature (i.e. the minimum slope in a logT vs. logi plot such 
as in Fig. 1) for several values of the inelastic parameter 7. As the collisions become more elastic, the decay rate approaches 
the theoretically predicted value of —2 for a homogeneous cooling state. Both axes are logarithmic, and the line is a power-law 
fit, with a slope of 0.3. 
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FIG. 3: The time at which the decay rate of the temperature reaches maximum, as a function of 7. The errors are rounding 
errors due to the finite sample rate. 
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FIG. 4: A flow diagram for 7 = 50 at t = 450. Each arrow represents the average velocity of four lattice sites, to improve 
readability. 

Figure 1 shows the decay of the average temperature as a function of time in our simulation, for three different 
values of the inelasticity parameter 7. The initial decay approximately follows the predicted —2 exponent, while for 
later times the temperature decays at a slower rate as the instabilities agitate the system[15, 16]. In the limit of 
low inelasticity, the maximal rate of decay more closely approaches Haff's predicted inverse-square behavior (Fig. 2). 
(Note, however, that the temperature will not decay at all in a completely elastic system.) In more inelastic systems, 
the hydrodynamic instabilities kick in sooner and compromise the homogeneity of the system. One could characterize 
the time it takes for the instabilities to emerge by the time it takes for the temperature to reach its fastest decay 
rate. Figure 3 shows that this onset time decreases with respect to the inelasticity parameter as a power-law with an 
exponent of —4/3. 

The first instability that is predicted to dominate the homogeneous solution is a hydrodynamic shearing mode: two 
bands of material moving in opposite directions. This has been seen in several molecular dynamic simulations [14, 17], 
Fig. 4 shows how it has developed in our system as two horizontal shear bands. Figure 5 shows that our system 
develops a clustering instability as well. Note that the single cluster takes on a compact shape, which is surprising 
given that there is no surface tension in our model. If clusters are supposed to grow by accretion, then one would 
expect a ramified structure. This puzzle becomes clearer as we consider a larger system after further evolution. 
Figure 6 shows that clusters first form in these compact shapes, but as time goes on they reach out to their neighbors, 
stretching into the more stringlike forms seen in simulation [4]. In hindsight, we are able to see this behavior in the 




FIG. 5: The density distribution for the 7 = 50 system at time t — 100. 

smaller system as well, where the single cluster interacts with itself through the periodic boundaries. 

Finally, to demonstrate that this clustering instability is the result of the inelastic parameter, we compare the width 
of the density distribution for inelastic systems with the distribution for the elastic case 7 = (Fig. 7). In the absence 
of inelasticity, the density distribution collapses to a delta function, indicating complete homogeneity. 

Our results show that the shearing and clustering instabilities, identified by Goldhirsch and Zanetti using a simplified 
version of the above equations, exist in the complete nonlinear granular hydrodynamic equations 1. Haff's cooling 
law is obeyed in the limit of small inelasticity, but in general the instabilities become relevant before the system 
has a chance to completely homogenize. The power-law dependence of the onset time for these instabilities on the 
inelasticity, and the —4/3 exponent in particular, are interesting; we have not found any reference to these in the 
literature. Also interesting is the way in which these clusters develop from compact structures into networks which 
span the system. It is not clear whether an individual cluster begins to stretch out because of the proximity of its 
neighbors, or because of effects due to its increasing surface size. One can imagine that the behavior of this system 
could change as we alter the total amount of mass in the system. 

We thank Professor Todd DuPont for his assistance in constructing our simulations, and Professors Heinrich Jaeger, 
Sidney Nagel, and Thomas Witten for helpful conversations. This work was supported by the Materials Research 
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FIG. 7: The width of the density distributions (i.e. p m ax ~ pmin) for several values of 7, including the elastic case 7 = 0. Notice 
that the density distribution is collapsing to a delta function in the elastic case, approaching complete homogeneity, while the 
inelastic systems show broadening density distributions due to the clustering instability. 
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